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Climatological  fields  of  mixed  layer  depth  ( MLD )  are  presented  over  the  Aegean,  Marmara,  Black  and  Azov  Seas. 
Monthly  fields  of  MLD  are  formed  by  historical  individual  temperature  and  salinity  profiles  from  combination 
of  various  data  sets  with  additional  quality  control  procedures  applied.  Various  definitions  that  are  based  solely 
on  temperature  (T)  or  those  that  include  the  impact  of  salinity  (S)  are  applied  to  investigate  the  robustness  in 
the  pattern  and  values  of  the  MLDs.  Interpolation  of  the  MLD  fields  to  a  0.25°  x  0.25°  regular  grid  over  the  region 
is  accomplished  using  a  combination  of  median  filter  and  ordinary  kriging.  Strong  seasonal  variability  is  noted 
in  all  regions.  Given  a  density-based  MLD  criterion  that  includes  both  T and  S,  deep  mixed  layers  (>200  m)  are 
noted  in  the  Aegean  Sea,  especially  eastern  part  of  the  region  during  winter  while  MLDs  are  generally  much 
shallower  (<60  m)  in  the  Black  Sea.  A  criterion  based  on  curvature  method,  which  determines  MLD  according 
to  first  maximum  of  curvature  of  T  may  fail  in  representing  deep  MLDs  in  the  Aegean  Sea  when  the  water 
column  is  well  mixed.  MLD  fields  obtained  from  all  definitions  are  found  to  be  strongly  correlated  to  each  other 
over  the  seasonal  cycle,  confirming  the  strong  seasonal  cycle.  While  the  curvature  method  gives  shallow  MLDs 
only  during  winter,  it  has  relatively  large  skill  in  comparison  to  the  density-based  MLD  criterion.  In  general, 
MLD  fields  suffer  from  lack  of  input  T  and  S  profiles  in  the  Marmara  and  Azov  Seas,  thus  they  may  not  be 
representative.  Monthly  MLD  fields  presented  in  this  paper  are  available  for  various  applications,  such  as  mixed 
layer  studies,  ocean  biology,  ocean  modeling  and  acoustic  propagation. 
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1.  Motivation 

There  have  been  various  studies  examining  major  features  of 
mixed  layer  depth  (MLD)  over  the  global  ocean  (e.g.,  Kara  et  al.,  2003; 
Lorbacher  et  al.,  2006).  There  are  also  some  regional  studies,  including 
the  North  Pacific  Ocean  (Ohno  et  al.,  2004;  Suga  et  al.,  2004),  the 
equatorial  Pacific  Ocean  (Lukas  and  Lindstrom,  1991),  the  western 
equatorial  Atlantic  (Pailler  et  al.,  1999),  the  South  China  Sea  (Aijun 
et  al.,  2006;  Qu  et  al.,  2007),  the  Mediterranean  Sea  (D'Ortenzio  et  al., 
2005)  and  the  Bering  Sea  (Miura  et  al.,  2002). 

Little  or  almost  no  attention  has  been  given  to  characteristics  of 
MLD  in  the  Aegean,  Marmara  and  Black  Seas.  As  part  of  the 
Mediterranean  Sea,  the  study  by  D'Ortenzio  et  al.  (2005)  included 
analysis  of  MLD  in  the  Aegean  Sea,  but  it  did  not  present  spatial  and 
temporal  variabilities  in  the  whole  region,  in  detail.  As  of  this  writing. 
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there  has  been  no  study  which  describes  seasonal  characteristics  of 
MLD  in  the  Aegean,  Marmara,  Black  and  Azov  Seas,  whose  bottom 
topography  is  quite  variable  (Fig.  1).  A  lack  of  the  understanding  of 
major  features  of  MLD  in  these  regions  has  motivated  this  study.  In 
addition,  the  knowledge  of  MLD  is  particularly  important  for  biological 
studies  in  the  Black  Sea,  where  the  autumn  and  spring  blooms 
usually  take  place  in  the  mixed  layer  (Oguz  et  al.,  2000),  and  oxygen 
distribution  is  of  particular  importance  for  biogeochemical  processes 
(Oguz  et  al.,  2001).  The  Marmara  Sea  is  also  of  particular  importance 
since  the  strong  density  gradient  layer  between  water  masses  may 
have  quite  large  influences  on  mixing  processes  (Besiktepe  et  al.,  1994). 

The  major  goal  is  to  present  spatial  and  temporal  variabilities  of 
MLD  on  the  climatological  time  scales  for  various  purposes,  such  as 
ocean  mixed  layer  studies,  ocean  biology,  near  surface  acoustic 
transmission  and  model-data  validations  of  MLD  from  an  ocean 
general  circulation  model  (OGCM).  In  studying  the  general  features  of 
climatological  monthly  means  of  MLD,  we  process  historical  tempera¬ 
ture  and  salinity  profiles  from  very  comprehensive,  up  to  date  and 
quality-controlled  data  sets,  extending  up  to  2006.  MLD  is  then 
determined  at  the  locations  where  individual  profiles  are  located.  This 
is  accomplished  with  various  definitions.  The  reason  for  applying 
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Fig.  1.  Water  depths  in  the  Aegean,  Marmara,  Black  and  Azov  Seas.  The  bottom  topography  is  constructed  from  the  1  min  resolution  data  from  General  Bathymetric  Chart  of  the 
Oceans  (GEBCO),  which  is  available  online  at  http://www.bodc.ac.uk/products/bodc_products/gebco/.  In  the  rest  of  the  paper,  we  use  this  domain  in  examining  MLD  feature  and  drop 
the  latitude  and  longitude  labels  from  figures  for  simplicity. 


different  MLD  definitions  based  on  T  only  or  including  the  impact  of  5 
through  density  is  to  examine  (1 )  how  robust  the  results  are,  (2)  how 
important  the  effect  of  salinity  is,  and  (3 )  what  relationship  might  exist 
among  the  MLD  fields. 

After  obtaining  the  irregularly-spaced  MLD  fields  at  individual  T 
and  5  locations  based  on  various  definitions,  fine  resolution 
(0.25° X 0.25°)  gridded  fields  are  produced  over  the  entire  region, 
involving  the  Aegean,  Marmara,  Black  and  Azov  Seas.  The  gridding  is 


done  based  on  a  method  for  interpolating  MLD  values  from  sample 
data  using  regionalized  variable  theory  in  which  the  prediction 
weights  are  derived  from  a  fitted  variogram  model,  i.e.,  ordinary 
kriging  (e.g.,  Diggle  and  Ribeiro,  2007). 

This  paper  is  organized  as  follows.  Section  2  describes  T  and  S 
profiles  in  forming  monthly  mean  MLD  fields  over  the  Aegean, 
Marmara,  Black  and  Azov  Seas.  Section  3  gives  various  MLD  definitions 
used  throughout  the  paper.  Section  4  provides  an  overview  of  spatially 


Fig.  2.  Locations  of  profiles  used  for  constructing  the  MLD  climatology  in  July:  (a)  temperature-only  profiles  and  (b)  temperature  and  salinity  profiles. 
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and  temporally  varying  MLD  fields  at  a  grid  resolution  of 
(0.25°  X 0.25°).  Section  5  discusses  differences  in  MLD  fields  obtained 
from  various  definitions.  Section  6  concludes  the  paper. 

2.  Data 

Observed  profiles  of  historical  temperature  and  profile  pairs  of 
temperature  and  salinity  from  three  sources  are  used  to  determine 
MLD:  These  are  as  follows:  (1)  The  World  Ocean  Database  2005 
(WOD05)  (Boyer  et  al.,  2006),  (2)  the  U.  S.  Navy's  Master  Oceano¬ 
graphic  Observation  Data  Set  (MOODS)  (Teague  et  al.,  1990),  and  (3) 
the  ARGO  float  data  (Gould  et  al.,  2004).  Currently,  no  temperature 
and  salinity  profiles  are  available  from  the  MOODS  archive  in  the 
Marmara  and  Azov  Seas.  Since  data  sources  for  MOODS  and  WOD05 
are  mostly  identical,  the  bulk  of  these  data  is  the  same.  In  addition,  the 
WOD05  has  less  stringent  quality  control  (QC)  procedures  than 
MOODS  and  potentially  more  data  sources.  For  these  reasons,  there  is 
about  a  6%  increase  in  data  volume  with  the  addition  of  WOD05  data 
relative  to  the  MOODS  data  set  alone. 

As  an  example,  locations  of  all  historical  profiles  are  marked  in 
Fig.  2  for  temperature-only  and  the  pairs  of  temperature  and  salinity. 
The  former  profiles  can  be  used  for  a  temperature-based  MLD  deter¬ 
mination,  and  the  latter  ones  allows  computation  of  density  and 
thereby  a  density-based  MLD.  The  pairs  of  temperature  and  salinity 
profiles  (958)  are  approximately  three  times  less  than  the  tempera¬ 
ture-only  profiles  (3188).  There  are  some  regions  where  either 
temperature  or  salinity  profiles  are  missing.  Thus,  there  are  fewer 
pairs  of  profiles. 

Table  1  lists  the  total  number  of  stations  over  three  basins, 
including  whole  Aegean,  Marmara  and  Black  Seas.  The  largest  number 
of  profiles  for  the  whole  region  occurs  in  May  (4765 -h  4429  =  9194 
total)  with  the  fewest  in  January  (1928-1-3137  =  5065  total).  The 
numbers  are  from  data  on  the  original  observation  depths  (not 
standard  depths)  before  the  quality  control  steps  described  in  Section 
4.  For  example,  the  numbers  of  profiles  listed  in  Table  1  with  regard  to 
Fig.  2  is  obtained  after  the  quality  control  steps.  In  particular,  the 
super-observation  step  eliminates  many  profiles  that  were  taken  at 
the  same  location  and  month  but  in  different  years.  The  table  entries 
labeled  as  high  resolution  are  profiles  that  do  not  have  any  vertical 
gaps  larger  than  20  m  in  the  upper  150  m  of  the  profile.  There  are 
relatively  less  data  at  high  vertical  resolution. 

Numbers  in  Table  1  are  also  given  for  the  eastern  and  western  parts 
of  the  Black  Sea,  separately.  In  addition  to  the  number  of  profiles  given 
in  the  table,  we  use  additional  14,145  profiles  in  the  Azov  Sea  for  the 


period  of  1913  through  2004.  The  month  with  the  greatest  number  of 
profiles  is  June  with  more  than  1800  and  the  month  with  the  fewest  is 
February  with  less  than  800.  Many  stations  have  very  few  vertical 
levels  and  will  not  be  used  to  compute  MLDs. 

We  should  emphasize  that  for  MLDs  that  require  only  temperature 
the  total  number  of  profiles  can  be  used.  Based  on  Table  1,  there  are 
enough  profiles  to  construct  a  climatology.  In  addition,  we  have  made 
considerable  efforts  to  minimize  the  effect  of  non-uniform  sampling, 
which  could  create  structure  that  is  not  physically  present,  using  the 
quality  control  procedures  described  in  Section  4. 

The  Marmara  and  Azov  Seas  are  two  regions  where  there  are  not 
many  profiles  to  compute  the  MLD  (Fig.  2).  This  poses  problems 
especially  when  a  density  criterion  is  used  for  computing  MLD  due  to 
relatively  less  number  of  salinity  profiles.  In  fact,  both  the  Marmara 
and  Azov  Seas  do  have  many  temperature  and  salinity  data  but  they 
are  mostly  near  the  sea  surface  down  to,  say  2  or  3  m,  such  as  those 
from  bottle  measurements.  MLD  is  not  computed  unless  there  are  at 
least  three  depth  levels  in  the  profile. 

While  the  profile  data  we  use  in  this  paper  are  quality-controlled, 
additional  errors  still  exist.  For  example  tens  of  profiles  in  the  Azov  Sea 
region  appear  to  have  their  latitude  and  longitude  locations  reversed 
putting  the  observations  over  land.  It  is  plausible  that  other  similar 
errors  occur  where  observations  from  other  oceans  are  mistakenly  put 
in  our  study  region.  Reversing  the  sign  of  the  longitude  for  example 
could  put  an  Atlantic  Ocean  profile  in  the  Black  Sea.  Many  of  these 
types  of  errors  have  been  caught  in  the  quality  control  procedures  of 
the  National  Oceanic  Data  Center  (NODC)  or  the  U.S.  Navy's 
Oceanographic  Office  but  the  state  of  the  art  of  quality  control 
procedures  is  not  perfect.  In  the  case  of  a  longitude  sign  error,  this  may 
only  be  caught  if  the  depth  of  the  profile  exceeded  the  bathymetric 
depth  for  the  recorded  location.  A  short  profile  in  the  Atlantic  Ocean 
could  then  be  misplaced  to  the  Black  Sea  by  passing  the  bathymetry 
check.  To  minimize  these  errors  in  the  resulting  MLD  fields,  outliers  are 
removed  statistically  as  will  be  described  in  Section  4,  in  detail. 

3.  MLD  methodologies 

Throughout  the  rest  of  the  text,  we  use  T,  S  and  (Jt  to  denote 
temperature,  salinity  and  density.  Note  that  (Jt  =  p  — 1000  kg  m“^, 
where  p  is  the  density  of  water.  The  density  profile  is  computed  using 
T  and  S  values  at  given  depths  based  on  the  standard  United  Nations 
Educational,  Scientific,  and  Cultural  Organization  (UNESCO)  equation 
of  state  with  no  pressure  dependence,  i.e.,  zero  pressure  (Millero  and 
Poisson,  1981).  The  inclusion  of  pressure  effects  increases  the  density 


Table  1 

Total  number  of  profiles  from  WOA05,  MOODS  and  ARGO  float  data  for  three  basins  (whole,  eastern  Black  Sea  and  western  Black  Sea)  by  data  source:  CTD  represents  profiles  with  T 
and  S  and  XBT  represents  profiles  with  T  only. 


Jan 

Feb 

Mar 

Apr 

May 

Jun 

Whole  region:  (22°E-4184°E,  34°N-47.3°N) 

CTD 

928 

2372 

3124 

2973 

4765 

2645 

XBT 

3137 

2633 

3594 

2822 

4429 

3631 

CTDhi 

167 

356 

555 

537 

792 

435 

XBThi 

1778 

1300 

1566 

1007 

1346 

1023 

Eastern  Black  Sea:  (34°E-4VE,  4VN-45°N) 

CTD 

464 

950 

1076 

1355 

2122 

1163 

XBT 

101 

128 

127 

143 

210 

222 

CTDhi 

64 

98 

85 

104 

214 

123 

XBThi 

13 

1 

29 

14 

11 

45 

Western  Black  Sea:  (28°E-34°E,  4VN-45°N) 

CTD 

130 

907 

1003 

1018 

1415 

702 

XBT 

1531 

861 

541 

187 

516 

586 

CTDhi 

9 

63 

74 

135 

205 

115 

XBThi 

1324 

623 

267 

9 

82 

73 

Jul 

Aug 

Sep 

Oct 

Nov 

Dec 

2068 

1901 

1502 

2325 

3235 

1902 

4040 

3649 

3787 

4162 

3440 

2550 

320 

305 

479 

455 

508 

342 

897 

850 

991 

763 

859 

487 

729 

540 

544 

895 

1185 

810 

462 

297 

295 

275 

198 

268 

85 

60 

130 

96 

124 

71 

30 

22 

17 

71 

24 

8 

652 

389 

344 

727 

1363 

613 

990 

742 

583 

369 

347 

462 

79 

65 

79 

102 

149 

121 

74 

77 

144 

53 

41 

29 

CTDhi  and  XBThi  represent  T and  S  and  T-only  profiles  respectively,  that  have  no  vertical  gaps  larger  than  20  m  in  the  upper  150  m  and  are  considered  “high  resolution.”  Note  that  the 
CTDs  and  XBTs  could  be  from  moorings  or  other  types  of  platforms.  In  the  text  they  are  refereed  to  as  T  and  S  profiles  and  T-only  profiles,  respectively. 
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Table  2 

Various  MLD  definitions  based  on  T  and  cij  through  T  and  5  used  in  this  paper. 


The  algorithms  for  the  first  three  ones  are  available  online  at  the  web  page  http:// 
www7320.nrlssc.navy.mil/nmld/nmld.html,  and  the  last  one  (i.e.,  MLD(L))  is  available 
at  http://www.ifm-geomar.de. 


gradient  sufficiently  rapidly  with  depth  that  it  produces  a  markedly 
shallower  MLD  that  is  strongly  inconsistent  with  the  MLDs  inferred 
from  the  corresponding  T  and  5  profiles  (e.g.,  Kara  et  al,  2000). 

One  of  the  goals  of  this  study  is  to  investigate  the  differences  in 
MLDs  calculated  based  on  T  only  and  T  and  5  in  Section  4.  Thus,  we 
determine  MLD  based  on  various  definitions  as  summarized  in  Table  2. 
MLD(crT)  refers  to  the  depth  at  the  base  of  an  isopycnal  layer  where 
the  density  has  changed  by  an  amount  of 


Aa^  =  a^{T  +  AT,S,P)  -a^{T,S,P),  (1) 

where  P=0,  from  the  density  at  a  reference  depth  of  3  m.  This 
definition  takes  full  account  of  density  changes  due  to  T  and  S 
variations  with  location  as  described  in  Kara  et  al.  (2000).  Unlike  the 
earlier  studies  using  a  constant  density  change  of  0.125  kg  m“^  (e.g., 
Suga  et  al.,  2004;  Ohno  et  al.,  2004),  the  MLD(ctt)  definition  involves  a 
variable  density  criterion  with  AT  value  of  0.8  °C.  While  this  specific  AT 
value  is  typically  suitable  for  mid-latitude  locations  (e.g.,  Kara  et  al., 
2000)  and  subject  to  tuning,  it  is  not  our  major  focus  here.  The  variable 
density  definition  is  appropriate  because  T  generally  compensates  for 
a  variable  density. 

The  MLD((Jc)  definition  is  the  depth  at  the  base  of  an  isopycnal 
layer,  where  the  Ot  value  differs  from  the  value  at  the  first  level  below 
3  m.  It  uses  a  constant  threshold  aj  value  equal  to  0.125  kgm“^.  On 
the  contrary,  MLD(ot)  is  estimated  by  searching  for  the  first  depth 
below  3  m  that  differs  from  a  variable  threshold  (Jt  value,  which  is 
computed  based  on  T  and  5  profiles.  Specifically,  the  searching 
threshold  is  equal  to  Oj  computed  at  zero  pressure,  the  reference  T 
minus  0.8°  and  the  reference  S.  The  threshold  is  therefore  variable  but 
corresponds  to  a  0.8°  temperature  change. 

In  the  case  of  MLD(r),  only  T  profiles  are  used.  Simply,  it  is  the 
depth  at  the  base  of  an  isothermal  layer,  where  the  Thas  changed  by  a 
fixed  amount  of  AT  =  0.8  °C  from  the  temperature  at  a  reference  depth 
of  3  m.  Such  definition  generally  represents  an  isothermal  layer  depth 
rather  than  MLD(o't).  The  difference  between  the  two  is  called  barrier 
layer  (e.g.,  Sprintall  and  Tomczak,  1992). 

Finally,  another  definition  (MLD(L))  as  introduced  by  Lorbacher 
et  al.  (2006)  is  applied.  The  algorithm  is  publicly  available  and  used 
here  as  it  is.  It  determines  MLD  by  finding  the  first  maximum  of 
curvature  of  T  profiles.  MLD(L)  has  the  advantage  that  the  estimated 
layer  depth  is  not  strongly  a  function  of  another  dynamical  quantities, 
such  as  the  surface  temperature.  Unlike  the  other  gradient-based 
methodologies,  the  methodology  suggested  by  Lorbacher  et  al.  (2006) 
does  not  strongly  depend  on  seasonal  and  regional  differences. 

The  threshold  methodology  (e.g.,  MLD(o't))  accommodates  the 
wide  variety  of  density  and  temperature  stratifications,  such  as  a 
subsurface  mixed  layer  underlying  a  surface  thermal  inversion; 
multiple  fossil  layers  beneath  the  surface  mixed  layer.  The  threshold 
methodology,  however,  is  limited  in  that  it  amounts  to  a  bulk 
computation.  That  will  return  the  layer  most  closely  associated  with 
a  seasonal  MLD.  The  curvature  methodology,  on  the  other  hand,  is  a 
measure  of  the  most  near  surface  curvature  peak  and  therefore  always 
returns  the  shallowest  measurable  MLD,  which  for  high-resolution 
profiles  could  be  the  diurnal  MLD  as  explained  by  Helber  et  al. 
(in  press). 


As  an  application  of  all  MLD  definitions  mentioned  above,  we 
present  T  and  S  and  calculated  density  profiles  at  43.1  °N,  29.8°E 
located  in  the  Black  Sea  on  February  of  1982  (Fig.  3).  Based  on  these 
profiles  MLD(cjt)  and  MLD((Jc)  are  almost  identical  with  values  of 
46  m  and  48  m.  The  temperature-only  criterion,  i.e.,  MLD  (I),  results  in 
deeper  value  of  71  m.  The  curvature-based  MLD(L)  definition  gives  a 
value  of  46  m  at  this  particular  location,  which  is  consistent  with  MLD 
((Jt)  and  MLD((Jc).  Note  that  while  the  MLD  of  46  m  is  at  the  top  of  the 
thermocline,  the  MLD  of  71  m  is  at  the  bottom  of  the  thermocline.  The 
top  of  the  thermocline  is  what  one  would  consider  the  end  of  the  truly 
well  mixed  layer.  At  this  particular  location,  there  is  an  indication  of 
the  barrier  layer  (25  m)  since  MLD((Jt)  is  different  than  MLD(r). 


4.  Monthly  MLD  climatology 

As  mentioned  in  Section  2,  there  can  be  some  outliers  in  the  T  and 
S  profiles  from  WOD05,  MOODS  and  ARGO  floats.  In  addition,  the 
historical  observations  represent  an  unknown  underlying  process  and 
contain  unavoidable  noise  when  attempting  to  construct  a  climatol¬ 
ogy.  The  noise  can  be  the  result  of  various  factors,  for  example,  inter¬ 
annual  variability,  systematic  and  random  errors,  and  inhomogeneous 
sampling.  Thus,  we  first  explain  how  the  errors  in  resulting  MLD  fields 
due  to  improper  profiles  are  reduced  before  presenting  spatial 
variations  of  monthly  climatologies. 

There  are  three  steps  to  reduce  outliers  in  the  MLD  fields.  We  only 
provide  a  brief  description  as  follows.  The  first  one  is  the  Z-score 
statistics,  as  described  by  Lanzante  (1996),  which  is  applied  to 
normalize  MLD  values  using  the  bi-weight  mean  and  standard 
deviation.  For  example,  as  evident  from  Table  1,  the  large  relative 
drop  in  numbers  for  the  high-resolution  profiles  indicates  that  many  of 
the  profiles  used  in  the  analysis  have  relatively  large  gaps  in  the  upper 
150  m.  If  large  vertical  gaps  produce  a  very  large  MLDs  in  a  few  profiles 
randomly,  then  the  Z-score  statistics  removes  these  outliers. 

The  second  step  for  reducing  outlier  includes  climatological  super¬ 
observations,  representing  the  median  value  and  location  of  MLDs 
from  all  years  for  a  given  month  that  are  within  a  1/12°  circle.  Since 


Fig.  3.  Vertical  profiles  of  T,  S  and  density  at  43.1  °N,  29.8°E  on  11  February  1982.  The  T 
and  S  profiles  are  from  the  WOD05  data  set,  and  density  was  computed  with  UNESCO 
equation  of  state  with  no  pressure  dependence,  as  explained  in  the  text.  Profiles  were 
taken  at  a  local  time  of  12  GMT. 
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the  gridded  resolution  is  1/4°  there  will  be  up  to  four  super¬ 
observations  for  two  adjacent  gridded  data  points.  A  set  of  super¬ 
observations  are  computed  by  iterating  through  all  the  original 
observations  and  replacing  them  with  super-observations.  The  final 
step  is  to  apply  a  median  filter  (e.g.,  Perreault  and  Hebert,  2007).  This 
is  done  for  replacing  its  value  with  the  median  of  the  surrounding 
observations  but  without  altering  the  location  of  the  observations. 
Further  details  about  each  one  of  these  steps  are  provided  in  Helber 
et  al.  (2009-this  issue). 


The  method  of  kriging  is  applied  to  obtain  MLD  fields  at  a 
resolution  of  0.25°  x  0.25°.  Kriging  is  a  set  of  linear  regression  routines 
which  minimize  estimation  variance  from  a  predefined  covariance 
model  (Higgle  and  Ribeiro,  2007).  The  methodology  is  based  on  the 
assumption  that  irregularly  placed  MLD  values  are  treated  regionally. 
The  regionalized  MLD  then  becomes  intermediate  between  a  truly 
random  variable  and  a  completely  deterministic  variable.  The  kriging 
methodology  takes  the  spatial  correlations  into  the  account  for  MLDs 
near  each  other;  MLDs  that  are  widely  separated  are  statistically 


Fig.  4.  Monthly  mean  MLD  climatology  constructed  at  a  0.25°  x  0.25°  grid  resolution.  The  MLD  methodology  uses  a  variable  density  criterion  as  described  in  the  text.  We  use  T and  S 
profiles  from  both  WOA05,  MOODS  and  ARGO  data  archives  in  computing  MLDs. 
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independent.  Details  of  the  kriging  procedure  applied  in  this  study  are 
provided  in  the  appendix.  Note  that  the  approach  that  we  have  been 
using  is  to  pre-process  the  data,  so  that  all  observations  have  the  same 
error  level.  That  is  the  purpose  of  the  super-observation  and  median 
filter  steps  as  described  earlier. 

We  now  introduce  the  seasonal  cycle  of  MLD(o't)  over  the  Aegean, 
Marmara,  Black  and  Azov  Seas  (Fig.  4).  A  MLD  value  of  >200  m  is  seen 
in  the  Aegean  Sea  during  winter,  but  it  becomes  as  shallow  as  25  m 
during  summer.  Seasonal  variability  is  also  evident  in  the  Black  Sea. 
Shallow  MLDs  in  the  Black  Sea  can  be  attributed  to  the  attenuation  of 
solar  radiation  during  summer  months  (Kara  et  al.,  2005).  There  are 
also  some  significant  variations  in  MLD  in  the  Marmara  Sea  but  much 
less  variation  in  the  Azov  Sea,  where  water  depth  is  already  very 
shallow.  As  mentioned  earlier,  the  lack  of  T  and  S  profiles  in  these  two 
regions  in  comparison  to  others  also  limits  the  accuracy  of  derived 
MLD  features. 

MLDs  are  generally  deeper  by  100-150  m  in  the  Aegean  Sea  than  in 
the  Black  Sea.  Deep  MLDs  of  >200  m  seen  in  the  western  part  of  the 
Aegean  Sea  do  not  exist  anywhere  in  Black  Sea  at  all.  In  fact,  the  entire 
northern  Levant  and  the  southern  Aegean  region  represent  the 
locations  of  intermediate  and  deep  water  formation  and  depending 
on  the  severity  of  winter  conditions,  MLDs  can  be  as  deep  as  400- 
500  m.  The  MLD  of  200-300  m  is  a  typical  case  during  the  years  with 
weak  convection.  This  structure  has  also  some  spatial  variability.  For 


example,  the  Rhodes  region  with  a  permanent  cyclonic  gyre  has 
always  deeper  mixed  layer  than  regions  located  further  east. 

Unlike  the  Mediterranean  Sea,  there  is  generally  strong  density 
stratification  in  the  Black  Sea,  typically  limiting  the  depth  of  MLD  and 
resulting  in  shallower  values.  MLD  rarely  exceeds  50  m  in  the  Black 
Sea  during  winter,  and  is  very  shallow  (<20  m)  during  spring  and 
summer.  Starting  from  November,  the  deepening  in  MLD  picks  up  in 
the  Black  Sea.  Previously,  Kara  et  al.  (2008)  demonstrated  that 
precipitation  contributes  to  the  spring-summer  buoyancy  and  the 
existence  of  shallow  mixed  layer  formation.  This  results  in  reduced 
turbulent  mixing  in  the  Black  Sea.  Note  also  that  MLDs  in  the  Azov  Sea 
are  slightly  deeper  during  winter  in  comparison  to  other  seasons.  A 
striking  feature  of  the  monthly  MLD  fields  is  the  existence  of  Caucasus 
and  Sevastopol  anticyclonic  eddies  in  the  Black  Sea. 


5.  Differences  among  MLD  definitions 

Similar  to  the  MLD(ctt)  definition  that  includes  a  variable  density 
criterion,  we  apply  other  methodologies  listed  in  Table  2  to  produce 
monthly  mean  fields.  This  is  done  for  comparisons  of  fields.  Our 
purpose  is  to  find  out  whether  or  not  there  are  significant  differences 
in  the  MLD  fields  depending  on  the  criterion  used.  Note  that  MLD(o't) 
and  MLD((Jc)  are  computed  based  on  profiles  of  land  S,  while  MLD(r) 
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Fig.  5.  Spatial  variability  of  climatological  MLDs  in  January  as  obtained  from  four  different  methodologies  listed  in  Table  2. 
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and  MLD(L)  are  computed  using  profiles  of  T  only.  Processing  of  the 
monthly  fields  is  same  for  all  MLD  definitions. 

5.1.  Climatological  fields 

Climatologies  of  monthly  MLD  fields  computed  based  on  different 
definitions  are  presented  in  Fig.  5  for  January.  There  are  no  large  dif¬ 
ferences  in  MLD  fields  obtained  by  MLD  (ctt)  versus  MLD  (ctc)  and  MLD(r). 
Major  differences  arise  when  using  MLD(L)  which  computes  the  MLD 
based  on  the  first  curvature  peak  of  near  surface  density  or  T  profiles.  A 
comparison  of  mean  fields  based  on  MLD((Tt)  versus  MLD((Jc)  indicates 
that  the  use  of  the  variable  density  criterion  rather  than  the  constant 
density  criterion  does  not  change  MLD  significantly,  except  in  the 
Aegean  Sea  (Fig.  5a, b).  However,  shoaling  of  MLD  does  also  occur  at 
some  parts  of  the  Black  Sea  when  using  MLD(oc)  as  opposed  to  MLD(ot). 

The  impact  of  salinity  in  determining  the  MLD  is  evident  from 
comparisons  of  fields  based  on  MLD(o't)  versus  MLD(T)  as  shown  in 
(Fig.  5a,c).  The  impact  of  salinity  stratification  is  to  make  MLDs 
shallower.  This  is  evident  in  all  regions.  The  MLD(L)  gives  much 
shallower  depths  in  comparison  to  those  obtained  from  others  (Fig.  5d). 
This  is  particularly  true  for  July  in  the  Aegean  Sea  as  well  (Fig.  6). 

Some  noise  in  monthly  fields  also  exists  in  July,  regardless  of  which 
MLD  definition  is  used.  While  there  are  many  T  and  S  profiles  in  this 
particular  month  (see  Fig.  2),  the  number  of  data  points  alone  does  not 
ensure  that  the  field  is  sampled  well  enough.  It  often  occurs  that  in  one 


year  there  will  be  a  large  number  of  observations  in  a  certain  location 
and  then  no  observations  for  many  years  in  that  same  location.  If  that 
one  year  was  anomalous  the  result  is  a  bull's  eye.  There  may  be  many 
observations  but  the  field  may  still  be  poorly  sampled.  That  is  a  problem 
that  generally  exists  in  historical  observations.  ARGO  is  one  observa¬ 
tional  platform  that  is  helping  even  out  the  sampling  in  this  way. 

In  general,  MLD(L)  values  are  still  relatively  shallow  in  July  (Fig.  6). 
The  algorithm  by  Lorbacher  et  al.  (2006)  looks  for  an  isothermal  layer 
at  the  surface.  It  determines  what  is  isothermal  based  on  the  variance 
within  the  profile.  In  the  winter,  when  the  entire  water  column  is  well 
mixed  and  the  variance  of  the  profile  is  uniform  throughout,  the 
algorithm  cannot  find  an  isothermal  layer  with  variance  less  than 
deeper  in  the  profile.  In  that  case,  the  algorithm  results  in  very  shallow 
values  and  even  <5  m.  In  addition,  it  is  noted  that  the  original 
algorithm  for  MLD(L)  does  not  impose  a  minimum  depth  value  to 
determine  the  MLD.  This  may  result  in  MLD  values  of  <  1  m.  For  this 
reason,  and  to  be  consistent  with  other  definitions  we  specified  a 
minimum  MLD  value  of  3  m  in  finding  MLD(L)  values.  Lorbacher  et  al. 
(2006)  also  introduced  a  definition  which  is  similar  to  what  we  use 
here,  but  density  profiles  are  input  to  the  algorithm.  We  did  not  make 
further  comparisons  against  the  curvature-based  MLD(L)  algorithm 
for  the  purpose  of  this  paper. 

Overall,  one  thing  to  note  from  MLD  fields  in  July  is  the  existence  of 
a  few  bull's  eyes  which  are  due  to  sampling  variation.  In  particular,  the 
bull's  eyes  are  in  the  original  data.  They  are  due  to  geophysical  errors. 


(a)  Jul:  MLD(crT)  (c)  Jul:  MLD(o-c) 


(c)  Jul:  MLD(T) 


(d)  Jul:  MLD(L) 


Fig.  6.  The  same  as  Fig.  5  but  for  July. 
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Fig.  7.  Time  series  of  monthly  mean  climatological  MLD  as  obtained  from  various 
definitions  at  43°N,  30°E. 


i.e.,  they  arise  due  to  regional  variability  that  is  not  properly  sampled. 
They  mainly  occur  when  observations  are  made  during  a  relatively 
rare  event  that  is  not  observed  again.  The  median  filter  may  enhance 
these  by  preserving  the  edge  that  forms  around  them.  On  the  other 
hand,  at  least  for  the  case  of  Black  Sea,  the  bull's  eyes  could  be 
representing  different  quasi-permanent  eddies  sampled  during  the 
field  surveys.  There  are  plenty  of  such  features  within  the  Black  Sea, 
and  some  of  them  are  smoothed  out  during  the  construction  of 
climatology  and  some  are  not  and  stay  there  resembling  bull's  eyes. 

To  provide  further  insight  about  differences  in  MLD  fields  obtained 
from  the  definitions  listed  in  Table  2,  time  series  of  monthly  mean  MLDs 
are  examined  at  two  locations.  We  randomly  chose  two  locations,  one  of 
which  is  located  at  43°N,  30°E  in  the  western  Black  Sea,  and  the  other  of 
which  is  located  at  39°N,  24°E  in  the  western  Aegean  Sea.  Variations  in 
MLD  are  clearly  evident  at  43°N,  30°E  over  the  seasonal  cycle  (Eig.  7). 
MLD(r)  gives  deeper  depths  than  others  during  winter,  and  MLD(L)  has 
the  shallowest  values  over  the  seasonal  cycle.  The  most  striking  feature 
of  MLDs  is  evident  at  39°N,  24°E  where  MLD(L)  is  significantly  different 
than  others  from  January  through  April  (Eig.  8).  In  fact,  depths  obtained 
from  the  MLD(L)  definition  are  again  shallow  by  ;^100  m.  This  is 
because  the  water  column  is  relatively  uniform  for  the  entire  depth  of 
the  profiles,  and  the  MLD(L)  algorithm  is  unable  to  distinguish  an 
isothermal  layer  from  the  variance  of  the  entire  profile  thereby  assuming 
that  the  isothermal  layer  is  shallow  or  non-existent. 

The  southeastern  Aegean  Sea  typically  has  very  shallow  values  when 
using  the  curvature-based  MLD(L)  criterion,  but  relatively  deep  ones 
when  using  MLD(o't)  and  MLD(o'c).  As  an  example  of  why  MLD(L)  gives 
very  different  values,  we  examine  T,  S  and  density  profiles  at  37.2°N, 
23.7°E  on  4  January  1990  (Eig.  9),  demonstrating  an  example  of  a  well 
mixed  wintertime  structure.  MLD(o't)  and  MLD(o'c)  have  values  of 
238  m  and  205  m,  respectively.  These  are  much  shallower  than  400  m, 
which  is  what  the  temperature-only  criterion  of  MLD(r)  gives.  The 
shallowest  value  is  34  m  obtained  from  the  MLD(L)  definition,  which 
appears  to  be  more  realistic.  This  is  because  the  profiles  indicate  a 
slightly  warmed  thin  surface  layer  (with  surface  T  of  15.3  °C  and  Ot  of 
28.9  kg  m  on  top  of  a  previously  generated  homogeneous  layer. 
These  two  water  masses  are  separated  by  the  temperature  and  density 


Fig.  9.  Vertical  profiles  of  T,  S  and  density  at  37.2°N,  23.7°E  on  4  January  1990.  The  T and 
S  profiles  are  from  the  WOD05  data  set,  and  density  was  computed  with  UNESCO 
equation  of  state  with  no  pressure  dependence,  as  explained  in  the  text.  Profiles  were 
taken  at  the  local  time  of  15  GMT. 


differences  smaller  than  the  criteria  set  in  the  analysis.  Thus,  clearly  it  is 
possible  to  obtain  misleading  information  about  MLD  features  in  the 
region,  depending  on  using  either  a  curvature  method  or  threshold 
method. 

Since  the  temperature  range  of  this  profile  is  <1.5  °C  over  a  depth  of 
600  m,  this  can  be  considered  a  very  well  mixed  ocean.  At  a  depth  of 
^20  m  there  is  a  density  inversion.  Due  to  penetrative  convection  the 
effective  mixed  layer  depth  may  be  much  deeper  than  the  profile  itself. 
Eor  these  reasons,  none  of  the  MLD  methods  is  particularly  accurate.  The 
MLD(L)  method,  however,  is  deceptively  shallow.  Since  the  MLD(L) 
measures  the  curvature  variability  of  the  profile  it  tends  to  be  shallow  in 
profiles  that  are  well  mixed  because  it  is  unable  to  distinguish  an 
isothermal  layer  from  the  relative  uniformity  over  the  whole  profile. 

5.2.  Statistical  comparisons  of  MLD 

Monthly  mean  values  obtained  from  MLD((Jc),  MLD(r)  and  MLD 
(L)  are  compared  with  those  from  MLD (ctt),  which  is  based  on  the 
density-based  criterion  over  the  Aegean,  Black,  Marmara  and  Azov 
Seas  over  the  seasonal  cycle.  Our  purpose  is  to  quantify  differences 
among  MLD  values.  Eirst,  monthly  mean  differences  are  evaluated. 
Differences  are  formed  at  each  grid  point,  and  then  spatial  maps  are 
obtained  (Eig.  10).  Differences  in  MLD  with  respect  to  MLD(o't)  are 
typically  small  in  all  months,  especially  in  summer.  As  expected,  large 
differences  arise  during  winter  when  MLD  is  deep.  The  most  obvious 
feature  of  the  differences  in  MLD  values  is  that  the  MLD(L)  gives  much 
shallower  values  (>50  m)  than  MLD(o't)  in  the  Aegean  Sea  in  January 
and  April  (Eig.  10a,b).  MLDs  from  all  algorithms  agree  with  those  from 
MLD((Jt)  quite  well  during  July  and  October  (Eig.  10c,d). 

MLD  fields  from  MLD((Jc),  MLD(T)  and  MLD(L)  are  further 
evaluated  based  on  more  comprehensive  statistical  metrics.  Since 
winter  MLDs  are  usually  much  deeper  than  summer  MLDs,  and  there 
can  be  geographical  differences  in  MLDs  (i.e.,  deep  ones  in  the  Aegean 
Sea  and  shallower  ones  in  the  Black  Sea),  we  will  also  use  a  non- 
dimensional  statistical  metric  for  fair  comparisons  of  fields.  This  is 
needed  since  there  will  be  a  skewed  distribution  for  MLD  values,  i.e., 
much  larger  deviations  above  the  mean  than  below  the  mean.  This  can 
depend  on  the  season  and  time  of  the  year. 
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Fig.  10.  Differences  in  MLDs  calculated  using  the  variable  density  criterion  (i.e.,  MLD(ar))  in  comparison  to  IVILD(r)  and  MLD(L).  Monthly  mean  differences  are  shown  in  January,  April,  July  and  October.  In  all  cases,  MLD(ar)  is  subtracted  from 
others.  Note  that  color  palettes  are  different. 
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The  statistical  metrics  used  for  comparing  MLD  time  series  are  mean 
error  (ME),  correlation  coefficient  {R)  and  non-dimensional  skill  score 
(SS).  Let  X/  (i  =  1,2,...,  n)  be  the  set  of  n  reference  MLD(o't)  values,  and 
let  Yj  (z  =  1,2,. . .,  n)  be  the  set  of  corresponding  values  of  MLDs  from  each 
one  of  other  algorithms,  i.e.,  MLD(crc),  MLD(r)  and  MLD(L).  Also  letX  ( Y) 
and  ax  (oy)  be  the  mean  and  standard  deviations  of  the  reference  (other 
methodology)  values,  respectively.  Following  Murphy  ( 1995)  and  Wilks 
(1995),  the  preceding  statistical  measures  can  be  expressed  as  follows: 

ME  =  Y  -  x;  (2) 

R=lf2iXi-^iYi-Y)/(a^ay),  (3) 

i  =  l 


SS  =  R^  - 

R-iay/ax) 

2 

\Y-X)/of 

^cond 

^uncond 

In  the  time  series  comparisons,  n  is  equal  to  12,  i.e.,  we  have  monthly 
mean  MLD  time  series  at  given  grid  point  over  the  seasonal  cycle 


(January  through  December).  ME  is  obtained  by  subtracting  values 
from  MLD((Jt)  from  those  from  each  one  of  other  methodologies.  It 
simply  represents  climatological  mean  difference  with  respect  to  MLD 
values  obtained  from  the  MLD((Tt)  definition.  R  value  is  a  measure  of 
the  degree  of  linear  association  between  the  time  series. 

The  SS  in  Eq.  (4)  is  the  fraction  of  variance  explained  by  two  non- 
dimensional  biases  (conditional  bias,  Bcond,  and  unconditional  bias, 
^uncond)  which  are  not  taken  into  account  in  the  R  formulation  (see 
Eq.  (3)).  In  brief,  Buncond  (also  called  systematic  bias)  is  a  non- 
dimensional  measure  of  the  difference  between  the  means  of  the  time 
series,  while  Bcond  is  a  measure  of  the  relative  amplitude  of  the 
variability  in  the  two  data  sets.  Note  that  R^  is  equal  to  SS  only  when 
Bcond  and  Buncond  are  zero.  These  two  biases  are  never  negative.  SS  is 
1.0  for  perfect  agreement,  and  is  negative  for  Bcond  +  Buncond>B^- 

Fig.  11  provides  spatial  fields  of  ME,  R  and  SS,  all  of  which  were 
calculated  over  the  seasonal  cycle.  In  comparison  to  MLD(crT), 
climatological  mean  MLD  bias  for  MLD(o'c)  and  MLD  (I)  is  small 
within  ±  10  m  over  most  of  the  region  (Fig.  11a).  The  only  exception  is 
the  mean  biases  from  the  MLD(L)  criterion,  giving  too  shallow  values 
in  the  Aegean  Sea.  On  the  contrary,  seasonal  cycle  of  MLD  from  all 


Fig.  11.  Statistical  evaluations  between  the  pairs  of  MLD(ar)  versus  MLD(ac),  MLD(ar)  versus  MLD(r)  and  MLD(ar)  versus  MLD(L)  over  the  seasonal  cycle:  (a)  climatological  annual 
mean  bias  in  meters,  (b)  linear  correlation  coefficient  and  (c)  non-dimensional  skill  score.  Statistical  values  were  computed  using  monthly  mean  time  series  at  each  grid  point  over 
the  seasonal  cycle. 
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definitions  agrees  with  each  other  well.  This  is  evident  from 
correlation  values  close  to  almost  1  all  over  the  region,  even  for  the 
MLD(L)  criterion  (Fig.  11b).  Thus,  MLD  seasonal  cycle  is  robust 
regardless  of  which  definition  is  preferred.  However,  such  agreement 
in  R  values  does  not  identify  the  differences  in  means  and  standard 
deviations  of  MLD  time  series  from  two  given  methodologies  as 
mentioned  above. 

Given  the  deep  (shallow)  MLD  value  during  winter( summer),  the 
non-dimensional  skill  score,  which  takes  conditional  and  uncondi¬ 
tional  biases  into  the  account,  provides  better  information  about  the 
relationship  between  MLD  pairs.  In  computing  skill  values  the 
logarithmic  transformation  is  applied  to  each  of  the  12  monthly  MLD 
fields  at  every  grid  point.  The  purpose  is  to  make  the  distribution  of 
deviations  more  symmetric  about  the  mean  since  the  distribution  is 
positively  skewed  (long  tail  for  larger  values),  i.e.,  deep  MLDs  during 
winter. 

The  perfect  skill  value  of  1  is  evident  between  MLD  (ctt)  versus 
MLD(cJc)  almost  everywhere  (Fig.  11c).  This  indicates  that  the  use 
of  constant  density  criterion  as  opposed  to  the  variable  one  yields 
almost  identical  results  in  this  particular  ocean  domain.  However, 
this  is  not  the  case  for  the  pair  of  MLD(o't)  versus  MLD(r)  in  the 
easternmost  part  of  the  Black  Sea  and  some  parts  of  the  Azov  Sea.  The 
MLD(L)  criterion  also  gives  results  very  close  to  the  MLD(o't)  criterion 
over  the  majority  of  the  domain  since  positive  skill  is  the  criterion  for 
a  minimal  level  of  acceptable  agreement.  However,  there  are  also 
negative  skill  values.  These  are  evident  in  the  northeastern  Aegean  Sea, 
all  Azov  Sea  and  a  few  small  locations  in  the  western  Black  Sea. 

6.  Summary  and  conclusions 

The  major  focus  of  this  paper  is  to  present  climatological  monthly 
mean  fields  of  MLD  in  the  Aegean,  Marmara,  Black  and  Azov  Seas. 
Temperature  and  salinity  profiles  compiled  from  various  data  sets  are 
used  to  produce  relatively  fine  resolution  (0.25°  x  0.25°)  MLD  values. 
Comprehensive  quality  control  procedures  are  applied  to  individual 
profiles  before  forming  the  monthly  climatologies.  Various  definitions, 
based  on  temperature-only,  density,  and  a  curvature  method;  were 
used  to  determine  MLD.  The  resulting  fields  were  also  compared  to 
each  other  to  determine  the  robustness  of  the  results. 

Overall,  all  definitions  reveal  the  existence  of  deep  (>200  m) 
mixed  layers  in  the  Aegean  Sea  during  winter.  On  the  contrary,  much 
shallower  MLDs  (<60  m)  are  noted  in  the  Black  Sea.  MLD  does  not 
change  significantly  during  summer  in  both  regions.  In  general,  while 
we  have  noted  a  few  differences  between  the  results  of  various  MLD 
definitions  over  the  whole  region,  the  annual  cycle  remains  typically 
same.  Not  much  seasonal  variability  is  seen  in  MLD  over  the  Marmara 
and  Azov  Seas  where  there  are  not  many  T  and  S  profiles,  greatly 
limiting  our  ability  to  obtain  reliable  fields  in  some  months.  Thus, 
more  accurate  determination  of  seasonal  MLD  variability  in  these 
regions  is  subject  to  future  studies.  In  particular,  the  Azov  Sea  is  not  an 
ideal  region  for  mixed  layer  studies  since  it  is  very  shallow,  and 
generally  <20  m  in  depth.  It  should  be  noted  that  in  summer  the  Azov 
Sea  is  strongly  stratified  due  to  insolation  which  has  a  strong  effect 
due  to  very  low  transparency  of  its  water  and  river  discharges.  Thus, 
the  mixed  layer  may  hit  the  bottom. 

In  this  paper,  we  do  not  make  any  concluding  remarks  about  which 
method  of  estimation  of  MLD  is  better  and  when  a  particular  one 
should  be  used.  This  requires  long  term  time  series  of  temperature  and 
salinity  at  some  particular  locations.  However,  it  is  demonstrated  that 
the  threshold  computation,  which  is  based  on  a  bulk  estimation,  can 
result  in  quite  different  MLDs  in  the  Aegean  Sea  in  comparison  to  the 
curvature  methodology.  The  latter  typically  returns  the  shallowest 
measurable  MLD.  Additionally,  our  results  reveal  that  both  methodol¬ 
ogies  generally  provide  similar  information  regarding  the  spatial  and 
temporal  variabilities  of  MLD  and  that  both  fail  under  certain 
circumstances.  For  example,  the  curvature  methodology  fails  when 


the  water  column  is  well  mixed  down  to  the  bottom  of  the  profile.  In 
this  case,  the  methodology  cannot  find  a  curvature  peak,  thereby 
returning  a  mixed  layer  value  at  the  shallowest  possible  point. 
Similarly  the  threshold  methodology  fails  when  there  is  a  weak 
shallow  stratification  with  temperature  or  density  deviations  smaller 
than  the  selected  threshold.  Evaluations  that  use  both  methodologies 
provide  a  more  robust  characterization  of  the  ocean. 
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Appendix  A 

Interpolation  of  the  irregular  MLD  values  to  0.25°  x  0.25°  grid  is 
performed  using  a  kriging  methodology.  The  kriging  takes  into 
account  the  presence  of  noise  and  constructs  a  linear  predictor 
based  on  the  covariance  structure  of  the  data.  The  software  used  for 
the  kriging  is  the  statistical  software  R  (http://www.r-project.org/) 
with  the  geostatistical  data  analysis  package  geoR  (http://www.leg. 
ufpr.br/geoR).  Details  about  the  package  can  be  found  in  Diggle  and 
Ribeiro  (2007). 

The  kriging  parameters  that  control  the  linear  predictor  and  the 
amount  of  smoothing  applied  to  the  output  field  are  the  range,  sill  and 
nugget.  The  predictor  model  is  spherical  with  the  parameters  and 
computed  by: 

(1)  the  range  value  of  (/)  =  2.5°, 

(2)  the  sill  value  for  nugget  =  0,  =  max(4.97(Jb  —  84.3, 10), 

(3)  nugget  =  1.4cr'o/(/>,  and 

(4)  cr'^  =  cr'^- nugget. 

The  zero  nugget  sill  value  of  is  computed  from  the  bi-weight 
standard  deviation  using  empirical  parameters  (4.97  and  —84.3) 
derived  from  a  linear  fit  with  the  sill  value  computed  from  a  likelihood 
fit  of  a  spherical  predictor  model  of  the  variogram.  Fitting  the  bi¬ 
weight  standard  deviation  through  the  empirical  parameters  is  done 
because  the  bi-weight  standard  deviation  is  a  more  robust  quantity 
than  the  variogram.  For  months  with  few  data  the  variogram  can  have 
unrealistic  values  while  the  bi-weight  standard  deviation  remains 
reasonable. 

The  kriging  field  is  further  smoothed  because  of  (3)  above  and  the 
limit  on  how  small  the  sill  can  be  (2).  The  range  value  is  based  on  the 
characteristic  length  scale  derived  empirically  from  variogram 
computations.  It  has  units  of  degrees.  Based  on  a  likelihood  fit  of  a 
spherical  predictor  model  to  the  observational  data,  the  mean  value  of 
the  sill  was  2.5°.  The  factor  of  1.4  in  (3)  is  chosen  to  represent  a 
characteristic  smoothing  factor  of  approximately  1.4°.  This  is  because 
cr'o/c/)  =  cr'o/ 2.5  =  nugget/ 1.4,  and  the  fact  that  2.5  is  the  character¬ 
istic  factor  based  on  the  variogram. 

Kriging,  as  we  are  using  it,  assumes  that  all  data  values  have  the 
same  error  level.  We  are  not  aware  of  a  kriging  procedure  that  would 
allow  for  an  error  specification  for  each  data  point.  Specifying  the 
error  level  for  each  data  point  could  help  weight  some  points  less  than 
others  and  thereby  even  out  data  errors  and  inhomogeneous 
sampling.  For  this  end,  an  inverse  solver  method  is  being  developed 
that  takes  into  account  variable  error  levels.  The  difficulty  is  in 
specifying  the  error.  It  may  be  that  specifying  the  error  is  not  realistic 
and  will  not  help.  One  of  our  other  future  goals  is  to  use  co-kriging 
methodology.  With  that,  we  could  use  different  combinations  of  MLD 
(ctt),  MLD (ctc),  MLD(r)  and  MLD(L)  to  get  a  more  robust  MLD  fields. 
One  could  also  use  the  temperature-only  estimates  to  supplement 
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the  T  and  5  estimates,  so  that  the  gaps  where  there  are  no  5  profiles 
can  be  filled. 
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